Beam propagation method

$$ i \frac{\partial \psi}{\partial z} = - \frac{1}{2 n_0 k_0} \nabla^2 \psi - \Delta n k_0 \psi $$$$ \frac{\partial \psi}{\partial z} = \frac{i}{2 n_0 k_0} \nabla^2 \psi + i \Delta n k_0 \psi $$$$ \frac{\partial \psi}{\partial z} = a \nabla^2 \psi + b \psi $$

where $ a = \frac{i}{2 n_0 k_0}$, $b = i \Delta n k_0$.

Discretization

$$ \nabla^2 \psi = \frac {\psi(x+h,y) + \psi(x-h,y) + \psi(x,y+h) + \psi(x,y-h) - 4 \psi(x,y)} {dx^2} $$

where $ h= dx$.

$$ \frac {\partial \psi}{\partial z} = \frac {\psi(z+dz;x,y) - \psi(z;x,y)} {dz}$$$$ \psi(z+dz;x,y) = \psi(z;x,y) + dz(a \nabla^2 \psi(z;x,y) + b \psi(z;x,y)) $$

Absorb Boundary

$$ \frac{\partial \psi}{\partial x}(z;-1,y) = \frac{\partial \psi}{\partial x}(z;1,y) = \frac{\partial \psi}{\partial x}(z;x,-1) = \frac{\partial \psi}{\partial x}(z;x,1) = 0 $$

CFL Condition

$$ \frac {a dz}{dx^2} \le \frac{1}{4} $$$$ dz \le \frac {dx^2}{4 a} $$$$ error \approx \frac{N}{3} (dz n_0 k_0)^3 $$$$ \frac {1}{n_0 k_0} (\frac{3}{N})^{-3} \le dz $$

Source Setting

$$ \frac {\partial \psi(z=z_0)} {\partial z} = f_0 $$$$ \psi(z=z_0) = f_1 $$



In [2]:

    
import numpy as np
import matplotlib.pyplot as plt
import math
%matplotlib inline



In [3]:

    
k0 = 2.0*math.pi/(633*10**-9)
n0 = 1.45
a = 1j/(2.0*n0*k0)
b = 0.0
L = 100.0*10*-6



In [7]:

    
size = 11*40  # size of the 2D grid
dx = L/size  # space step
Z = 10.0*10**-3
n = 10000
dz = Z/n;



In [8]:

    
psi = np.zeros([size, size],dtype=complex)
x = np.linspace(-L/2,L/2,size)
y = np.linspace(-L/2,L/2,size)



In [9]:

    
for i in range(size):
    for j in range(size):
        psi[i,j] = np.exp(-(x[i]**2 + y[j]**2)/(L/10)**2)
plt.imshow(abs(psi),cmap=plt.cm.hot);



In [10]:

    
def laplacian(Z):
    Ztop = Z[0:-2,1:-1]
    Zleft = Z[1:-1,0:-2]
    Zbottom = Z[2:,1:-1]
    Zright = Z[1:-1,2:]
    Zcenter = Z[1:-1,1:-1]
    return (Ztop + Zleft + Zbottom + Zright - 4 * Zcenter) / dx**2



In [13]:

    
# We simulate the PDE with the finite difference method.
for i in range(n):
    deltapsi = laplacian(psi)
    psic = psi[1:-1,1:-1]
    psi[1:-1,1:-1] = psic + dz * (a * deltapsi + b*psic)
    
    psi[0,:] = psi[1,:]
    psi[-1,:] = psi[-2,:]
    psi[:,0] = psi[:,1]
    psi[:,-1] = psi[:,-2]
    
    if i%(n/20) == 0:
        print 'Finished = {num}%'.format(num=100.0*i/n)
plt.imshow(abs(psi),cmap=plt.cm.hot);









    



Finished = 0.0%
Finished = 5.0%
Finished = 10.0%
Finished = 15.0%
Finished = 20.0%
Finished = 25.0%
Finished = 30.0%
Finished = 35.0%
Finished = 40.0%
Finished = 45.0%
Finished = 50.0%
Finished = 55.0%
Finished = 60.0%
Finished = 65.0%
Finished = 70.0%
Finished = 75.0%
Finished = 80.0%
Finished = 85.0%
Finished = 90.0%
Finished = 95.0%



In [203]:

    
plt.imshow(abs(psi),cmap=plt.cm.hot);



In [19]:

    
import numpy as np
import matplotlib.pyplot as plt
import math
%matplotlib inline

k0 = 2.0*math.pi/(633*10**-9)
n0 = 1.45; dn0 = 7.0*10**-4
a = 1j/(2.0*n0*k0)
b = 1j*k0

N = 10
A = 50.0*10**-6
L = A*N
La = 11.0*10**-6
Lb = 4.0*10**-6


size = N*50  # size of the 2D grid
dx = L/size  # space step
Z = 1.0*10**-3
n = 2000
dz = Z/n;

psi = np.zeros([size, size],dtype=complex)
dn = np.zeros([size,size])
x = np.linspace(0,L,size)
y = np.linspace(0,L,size)

for i in range(size):
    for j in range(size):
        psi[i,j] = np.exp(-((x[i]-A*N/2.0-A/4)**2/(1.5*A)**2 + (y[j]-A*N/2.0-A/4)**2/(1.5*A)**2))
        x0 = x[i]%A; y0 = y[j]%A
        dn[i,j] = dn0*(np.exp(-((x0-A/4)**2/(La/2)**2 + (y0-A/4)**2/(Lb/2)**2))\
                       + np.exp(-((x0-A*3/4)**2/(La/2)**2 + (y0-A/4)**2/(Lb/2)**2))\
                       + np.exp(-((x0-A/4)**2/(La/2)**2 + (y0-A*3/4)**2/(Lb/2)**2)))

psi = psi/np.sum(np.abs(psi))
plt.imshow(dn,cmap=plt.cm.hot)
plt.show()

def laplacian(Z):
    Ztop = Z[0:-2,1:-1]
    Zleft = Z[1:-1,0:-2]
    Zbottom = Z[2:,1:-1]
    Zright = Z[1:-1,2:]
    Zcenter = Z[1:-1,1:-1]
    return (Ztop + Zleft + Zbottom + Zright - 4 * Zcenter) / dx**2

for i in range(n):
    deltapsi = laplacian(psi)
    psic = psi[1:-1,1:-1]
    psi[1:-1,1:-1] = psic + dz * (a * deltapsi + b*dn[1:-1,1:-1]*psic)
    
    psi[0,:] = psi[1,:]
    psi[-1,:] = psi[-2,:]
    psi[:,0] = psi[:,1]
    psi[:,-1] = psi[:,-2]
    
    if i%(n/10) == 0:
         print 'Finished = {num}%, total field = {field} '.format(num=100.0*i/n,field=np.sum(np.abs(psi)))
plt.imshow(abs(psi),cmap=plt.cm.hot);









    












    



Finished = 0.0%, total field = 1.00000019091 
Finished = 10.0%, total field = 0.999305171709 
Finished = 20.0%, total field = 0.996430620806 
Finished = 30.0%, total field = 0.991380838707 
Finished = 40.0%, total field = 0.984942086262 
Finished = 50.0%, total field = 0.977810212173 
Finished = 60.0%, total field = 0.96740293406 
Finished = 70.0%, total field = 0.95257216782 
Finished = 80.0%, total field = 0.93521271287 
Finished = 90.0%, total field = 0.917408535053



In [18]:

    
import numpy as np
import matplotlib.pyplot as plt
import math
%matplotlib inline

k0 = 2.0*math.pi/(633*10**-9)
n0 = 1.45; dn0 = 7.0*10**-4
a = 1j/(2.0*n0*k0)
b = 1j*k0

N = 10
A = 50.0*10**-6
L = A*N
La = 11.0*10**-6
Lb = 4.0*10**-6


size = N*50  # size of the 2D grid
dx = L/size  # space step
Z = 1.0*10**-3
n = 2000
dz = Z/n;

psi = np.zeros([size, size],dtype=complex)
dn = np.zeros([size,size])
x = np.linspace(0,L,size)
y = np.linspace(0,L,size)

for i in range(size):
    for j in range(size):
        psi[i,j] = np.exp(-((x[i]-A*N/2.0-A/4)**2/(1.5*A)**2 + (y[j]-A*N/2.0-A/4)**2/(1.5*A)**2))
        x0 = x[i]%A; y0 = y[j]%A
        dn[i,j] = dn0*(np.exp(-((x0-A/4)**2/(La/2)**2 + (y0-A/4)**2/(Lb/2)**2))\
                       + np.exp(-((x0-A*3/4)**2/(La/2)**2 + (y0-A/4)**2/(Lb/2)**2))\
                       + np.exp(-((x0-A/4)**2/(La/2)**2 + (y0-A*3/4)**2/(Lb/2)**2)))

psi = psi/np.sum(np.abs(psi))
plt.imshow(dn,cmap=plt.cm.hot)
plt.show()

def laplacian(Z):
    Ztop = Z[0:-2,1:-1]
    Zleft = Z[1:-1,0:-2]
    Zbottom = Z[2:,1:-1]
    Zright = Z[1:-1,2:]
    Zcenter = Z[1:-1,1:-1]
    return (Ztop + Zleft + Zbottom + Zright - 4 * Zcenter) / dx**2

for i in range(n):
    deltapsi = laplacian(psi)
    psic = psi[1:-1,1:-1]
    psi[1:-1,1:-1] = psic + dz * (a * deltapsi + b*dn[1:-1,1:-1]*psic)
    
    psi[0,:] = psi[1,:]
    psi[-1,:] = psi[-2,:]
    psi[:,0] = psi[:,1]
    psi[:,-1] = psi[:,-2]
    
    if i%(n/10) == 0:
         print 'Finished = {num}%, total field = {field} '.format(num=100.0*i/n,field=np.sum(np.abs(psi)))
plt.imshow(abs(psi),cmap=plt.cm.hot);









    












    



Finished = 0.0%, total field = 1.00000019091 
Finished = 10.0%, total field = 0.999305171709 
Finished = 20.0%, total field = 0.996430620806 
Finished = 30.0%, total field = 0.991380838707 
Finished = 40.0%, total field = 0.984942086262 
Finished = 50.0%, total field = 0.977810212173 
Finished = 60.0%, total field = 0.96740293406 
Finished = 70.0%, total field = 0.95257216782 
Finished = 80.0%, total field = 0.93521271287 
Finished = 90.0%, total field = 0.917408535053



In [ ]: